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Abstract. In this paper we study the distribution function P(u a ) of the estimators 
u a ~ T _1 J o w(i)B 2 di, which optimise the least-squares fitting of the diffusion 
coefficient Df of a single d-dimcnsional Brownian trajectory B t . We pursue here 
the optimisation further by considering a family of weight functions of the form 
u)(t) = (to + t)~ Q , where to is a time lag and a is an arbitrary real number, and seeking 
such values of a for which the estimators most efficiently filter out the fluctuations. We 
calculate P(u a ) exactly for arbitrary a and arbitrary spatial dimension d, and show 
that only for a = 2 the distribution P(u a ) converges, as e = to/T —> 0, to the Dirac 
delta-function centered at the ensemble average value of the estimator. This allows us 
to conclude that only the estimators with a = 2 possess an ergodic property, so that the 
ensemble averaged diffusion coefficient can be obtained with any necessary precision 
from a single trajectory data, but at the expense of a progressively higher experimental 
resolution. For any a/2 the distribution attains, as e — > 0, a certain limiting form 
with a finite variance, which signifies that such estimators are not ergodic. 
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1. Introduction 



Single particle tracking (SPT) is an increasingly used method of analysis in biological 
and sot matter systems where the trajectories of individual particles can be optically 
observed. Recent advances in image processing, data storage and microscopy have led 
to an increasing number of papers, notably in biophysics, on single particle tracking 
in biological settings such as the cellular interior and the cell membrane. However 
the basic method of SPT owes its origin to the work of Perrin on Brownian motion 
[I], where optical observation is used to generate the time series for the position of an 
individual particle trajectory B t in a medium (see, e.g., Refs. [2 [3]). Complemented 
by the appropriate theoretical analysis, the information drawn from a single, or a 
finite number of trajectories, provides insight into the underlying physical mechanisms 
governing the transport properties of the particles. Via the analysis of the stochastic 
processes manifested in single particle trajectories, SPT is routinely used for the 
microscopic characterisation of the thermodynamic rheological properties of complex 
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media [I] , and also to identify non-equilibrium biological effects, for example the motion 
of biomolecular motors [5]. In biological cells and complex fluids, SPT methods have 
become instrumental in demonstrating deviations from normal Brownian motion of 
passively moving particles (see, e.g., Refs. [6l 171 151 l9l [TO] ) . The method is thus potentially 
a powerful tool to probe physical and biological processes at the level of a single molecule 

The reliability of the information drawn from SPT analysis, obtained at high 
temporal and spatial resolution but at expense of statistical sample size is not 
always clear. Time averaged quantities associated with a given trajectory may be 
subject to large trajectory-to-trajectory fluctuations. For a wide class of anomalous 
diffusions, described by continuous-time random walks, time-averages of certain 
particle's observables are, by their very nature, themselves random variables distinct 
from their ensemble averages [T21 H3]- For example, the square displacement time- 
averaged along a given trajectory differs from the ensemble averaged mean squared 
displacement |13[ [HI US]- By analyzing time-averaged displacements of a particular 
trajectory realization, subdiffusive motion can actually look normal, although with 
strongly differing diffusion coefficients from one trajectory to another [131 El and 
show pronounced ageing effects [T6j . 

Standard Brownian motion is a much simpler random process than anomalous 
diffusion, however the analysis of its trajectories is far from being as straightforward as 
one might think, and all the above mentioned troublesome problems exist for Brownian 
motion as well. Even in bounded systems, substantial manifestations of trajectory-to- 
trajectory fluctuations in first passage time phenomena have been revealed [T7J [18]. If 
one is interested in determining the diffusion coefficient of a given particle, standard 
fitting procedures applied to finite albeit very long trajectories unavoidably lead to 
fluctuating estimates Df of the diffusion coefficient, which might be very different from 
the true ensemble average value D, defined in a standard fashion as 

D= V , 1-1 

2dt ' V 1 

where the symbol E{. . .} denotes the ensemble average and d is the spatial dimension. 

To get a rough idea of how basic estimators for diffusion constants can fluctuate, 
consider a simple-minded, rough estimate oi Df, defining it as the slope of the line 
connecting the starting and the end-points B t of a given trajectory, i.e., like D in 
Eq. ( 11. ID but without averaging, that is Df = Bf/2dt. A single trajectory diffusion 
coefficient Df so defined is a random variable whose probability density function (pdf) 
P(Df) is the so-called chi-squared distribution with d degrees of freedom: 

P(Df) = -4r- ( Df~' exp L* ■ £l) , (1.2) 
v /; r(d/2) \2D / f V 2 D J 

where T(-) is the Gamma-function. The pdf in the latter equation diverges as Df — > 
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for d — 1, P{Df = 0) is constant for d = 2, and only for d > 2 the distribution has 
a bell-shaped form with the most probable value D*j = (1 — 2/d)D. This means that, 
e.g., for d = 3, it is most likely that extracting Dj from a single Brownian trajectory 
using this simple-minded approach we will obtain just one third of the true diffusion 
coefficient D. 

As a matter of fact, even taking advantage of more sophisticated fitting procedures, 
variations by orders of magnitude have been observed in SPT measurements of the 
diffusion coefficient for diffusion of the LacI repressor protein along elongated DNA 
[T9] . in the plasma membrane [3] or for diffusion of a single protein in the cytoplasm 
and nucleoplasm of mammalian cells [20]. Such a broad dispersion of the values of 
the diffusion coefficient extracted from SPT measurements raises important questions 
about the optimal methodology, much more robust than the rough estimate mentioned 
above, that should be used to determine the true ensemble average value of D from 
just a single trajectory. Clearly, it is highly desirable to have a reliable estimator even 
for the hypothetical pure cases, such as, e.g., unconstrained standard Brownian motion. 
A reliable estimator must possess an ergodic property so that its most probable value 
should converge to the ensemble average one and the variance should vanish as the 
observation time increases. This is often not the case and moreover, ergodicity of a 
given estimator is not known a priori and has to be tested for each particular form of 
the estimator. On the other hand, knowledge of the distribution of such an estimator 
could provide a useful gauge to identify effects of the medium complexity as opposed 
to variations in the underlying thermal noise driving microscopic diffusion. Recently, 
much effort has been invested in the analysis of this challenging problem and several 
important results have been obtained for the estimators based on the time-averaged 
mean-square displacement [221 E31 121] , mean maximal excursion [25] or the maximum 
likelihood approximation and its ramifications [261 [271 EHl EHl 130] . 

Let us define the dimensionless estimator of the diffusion coefficient as u = Df/D. 
In this paper, following our succinct presentation in [31], we focus on a family of least- 
squareijj estimators u a given by 

u a = ^£co(t)B 2 t dt, (1.3) 

where u(t) is the weight function of the form 

a being a tunable exponent, (positive or negative), to - a lag time and A a - the 
normalisation constant, appropriately chosen such that E{-u a } = 1. Therefore 

| This term will be made clear in Section [5] 
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Note that such a normalisation permits a direct comparison of the effectiveness of 
estimators corresponding to different values of a. Our goal here is to find such a for 
which u a is ergodic, namely, so that the single trajectory diffusion coefficient Df^-D 
(or u a — )• 1) as e = to/T — > 0. 

It should be emphasised that, as a matter of fact, we are dealing here with 
two consecutive optimisation schemes: first, the estimators in Eq. ( II. 3p are already 
the results of an optimisation of the least-squares fitting procedure for the diffusion 
coefficient Df of a single Brownian trajectory and second, an optimisation is performed 
for the weight function u(t) chosen from the family of functions in Eq. (jl.4p . 

This paper is outlined as follows: We start in Section [2]with a physical interpretation 
of the estimators given in Eq. (jl.3p and show that these stem out of the minimisation 
procedure of certain least-squares functionals of the square displacement B^. Next, in 
Section[3]we introduce basic notations and the definitions of the characteristic properties 
we are going to study. Section H] is devoted to the evaluation of the moment-generating 
function of the estimators. In Section [5] we present explicit results for the variance of 
the least-squares estimators for a ^ 2 for an infinite observation time, the variance 
for the case a = 2 for arbitrary observation time. We also discuss the optimisation of 
the variance of the least-squares estimators in the case of continuous-time and -space 
Brownian trajectories recorded at discrete time moments. Next, in Section [6] we obtain 
the asymptotic behaviour of the distribution P(u a ) in arbitrary spatial dimension. 
Further on, Section [7] discusses the shape of the full distribution P{u a ) and the location 
of the most probable values of the estimators for three and two-dimensional systems. 
We also study the distribution of the variable oj a = v$ /(u^ + «1 2 '), with Ua and u«' 
two independent estimates. Finally, in Section [8] we conclude with a brief summary of 
our results. 

2. Physical interpretation of the estimators u a . 

Before we proceed, it might be instructive to understand where do the functionals in 
Eq. (11.3P stem from and what physical interpretation may they have. Consider a given d- 
dimensional trajectory B< with t G [0, T] and try to calculate the diffusion coefficient Df 
of this given trajectory using the least-squares approximation for the whole trajectory. 
To this purpose, one writes the following least-squares functional 



and seeks to minimise it with respect to the value of Df, considered as a variational 
parameter. Note that Eq. (12. ip is a bit more general compared to the usually used 




(2.1) 
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least-squares functionals. A novel feature here is that in Eq. (12. ip we have introduced 
a weight function u(t) which, depending on whether it is a decreasing or an increasing 
function of t, will sensitive to the short time or the long time behavior of the trajectory 
Bj, respectively. 

Furthermore, setting the functional derivative OF/dDf equal to zero, we find that 
the minimum of F is obtained for 



Next, identifying the denominator with 1/A a in Eq. ( II. 5p . we conclude that u a in 
Eq. (jl.3p minimises the least-squares functionals with a weight function u(t) = (t +t)~ a . 

It is interesting to note that with this weight function, the functional ( II. 3p 
interpolates two well known estimators for the diffusion constant: In the case of a = — 1, 
the estimator u a corresponds to a usual unweighted least-squares estimate (LSE) of the 
time-averaged squared displacement [HI [201 [2T] . The case a = 1 arises in a conceptually 
different fitting procedure in which the conditional probability of observing the whole 
trajectory B t is maximised, subject to the constraint that it is drawn from a Brownian 
process with the mean-square displacement 2dDt, Eq. (II. ip . This is the so-called 
maximum likelihood estimate (MLE) which takes the value of Df that maximises the 
likelihood of B t , defined as: 

L = f[^D f tr^e W (-^-) , (2.3) 
t=o \ / / 

where the trajectory B t is appropriately discretized. Differentiating the logarithm of 
L with respect to Df and setting dlnL/dDf = 0, one finds the maximum likelihood 
estimate (see, e.g., Refs.[26 | l29 | l30]) of Df, which upon a proper normalisation is defined 
by Eq. flT3D with a = 1. 

3. Basic notations and definitions 

The fundamental characteristic property we will focus on is the moment-generating 
function $((r) of the random variable in Eq. f l 1 . 3 P : 

$(<t) = E {exp (-au a )} , (3.1) 

where a is a parameter. 

It is important to realise that for standard Brownian motion the generating function 
of the original <i-dimensional problem decomposes into a product of the generating 
function of the components, since the squared distance from the origin of a given 
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realisation of a d- dimensional Brownian motion at time t decomposes into the sum 



B? = $> t 2 «, (3-2) 



B t [i) being realisations of trajectories of independent ID Brownian motions (for each 
spatial direction). Thus, the moment-generating function $(<t) factorizes 

$(<r) = G d (a), (3.3) 

where 

G(a) = E jexp {-^ £ w (r) B^) dr^j } . (3.4) 

In what follows we will thus skip the argument (i) and denote as B t a given trajectory 
of a one-dimensional Brownian motion with an ensemble average diffusion coefficient D. 

The knowledge of <&(er) will allow us to calculate directly, by merely differentiating 
$(cr), the variance Vai(u a ) of the distribution function P(u a ) and to infer the asymptotic 
behaviour of the distribution function. The complete distribution P(u a ) will be obtained 
by inverting the Laplace transform in Eq. ( 13. ip with respect to the parameter a, namely 



P{u a ) = — / da exp (au a ) $(cx) , (3.5) 

where 7 is a real number chosen in such a way that the contour path of integration is 
in the region of convergence of $(cr). The explicit results for the variance and for the 
distribution P(u a ) will be presented in the next sections. 

Further on, to highlight the role of the trajectory-to-trajectory fluctuations, we will 
consider the probability density function P(u a ) of the random variable 

Ua = (1) ~ (2) ' ( 3 - 6 ) 

where Ua and are two identical independent random variables with the same 
distribution P{u a ). The distribution P(u a ), introduced recently in Ref.fT?] [T8] (see 
also Refs.(32j |33j [M]) is a robust measure of the effective broadness of P(u a ), and 
probes the likelihood that the diffusion coefficients drawn from two different trajectories 
are equal to each other. This characteristic property can be readily obtained from the 
identity [55] 
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Figure 1. (Color online) The distribution P(cj) in Eq. ([3~8]) for d = 1, 2 and 3. 



l-w, 



du a u a P(u a ) P 



0J r 



-u, 



(3.7) 



Hence, P(uj a ) is known once we know P(u a ). To illustrate this statement, let us return to 
the simple-minded estimate of Df mentioned in the Introduction and the corresponding 
pdf given by Eq. (11.21) . In this case, P(u) can be obtained explicitl}|§|, 



p ( w ' = f^p) w ' !/2 " 1(1 - w) ' i/2 "' (3 - 8) 

A striking feature of the beta-distribution in Eq. (13. 8p is that its very shape depends 
on the spatial dimension d (see Fig. [T|). For d = 1, P(oj) is bimodal with a ?7-like 
shape, the most probable values being and 1 and, remarkably, u = 1/2 being the least 
probable value. Therefore, if we take two Id Brownian trajectories, most likely we will 
obtain two very different estimators of the diffusion coefficient by this method, both 
having little to do with the true ensemble average value D. It is unlikely that we will 
get two equal values. Further on, for d = 2, P(u) = 1, which signifies that any relation 
between estimates Df drawn from two different trajectories is equally probable. Only 
for d = 3 the pdf P(oo) is unimodal with a maximum at u = 1/2. But even in this 
case it is broad (a part of a circular arc), revealing important trajectory-to-trajectory 
fluctuations. Clearly, such a simple-minded estimate is not plausible and one has to 
resort to more reliable estimators. Below we will consider the forms of P(oo a ) for the 
family of weighted least-squares estimators. 

§ We drop the subscript a as meaningless in this case. 
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4. The moment-generating function of the estimators. 

Note that the Laplace transforms of quadratic functionals of Brownian motion (and 
other Gaussian processes), as the one in Eq. (13. ip . have attracted a great deal of 
interest over the last decades following the pioneering work by Cameron and Martin [35] . 
Numerous results have been obtained both in the general setting of abstract Gaussian 
spaces and in various specific models generalising the original approach for Brownian 
motion due to Cameron and Martin (see, e.g., Refs. [36 ] I3T ] l38 | 139] for contributions and 
references therein). An alternative approach is based on the path integrals formulation 
for quantum mechanics, which represents a powerful tool to analyse these problems 
using methods more familiar to physicists [IUJ ST]: here, the problem appears as the 
computation of the partition function of a quantum-harmonic oscillator with time 
dependent frequency. Various quadratic functionals of Brownian motion have been 
studied in details by physicists [H] with a variety of methods. They arise in a plethora 
of physical contexts, from polymers in elongational flows [13] to Casimir/van der Waals 
interactions and general fluctuation induced interactions [MJ S3 SS, S3 SH] where, 
in the language of the harmonic oscillator, both the frequency and mass depend on 
time. Quadratic functionals of Brownian motion also arise in the theory of electrolytes 
when one computes the one-loop or fluctuation corrections to the mean field Poisson- 
Boltzmann theory [4*9| |50| |5T| [52] . Finally we mention that functionals of Brownian 
motion also turn out to have applications in computer science [53] . 

In order to calculate G(a) in Eq. ( 13 .411 we resort to the path integrals formulation. 
Following Refs. [291 EH] , we introduce an auxiliary functional 



*(*,*)= EN exp(-^ / uj{T)BldT)\ (4.1) 



where the expectation is for a Brownian motion starting at x at time t. In terms of 
ty(x,t) the moment-generating function is determined by noting that G(a) = \l/(0,0). 

Further on, we derive a Feynman-Kac type formula for ^f?(x, t) considering how the 
functional in Eq. (14. TJ evolves in the time interval (t, t + dt). During this interval the 
Brownian motion moves from x to x + dB t , where dB t is an infinitesimal Brownian 
increment such that EdB{dB t } = and KdB{dBf} = 2Ddt, where K^b denotes now 
averaging with respect to the increment dB t . For such an evolution we have, to linear 
order in dt 

*(*, t) = E dB { (l - ^f^x^ E£J* {exp jf W (r) B\ dr) } } 

= E dB S^(x + dB u t + dt) (l - aAc £® x 2 dt) | • (4.2) 

Expanding the right-hand-side of the latter equation to second order in dB t and to linear 
order in dt, we eventually find, after averaging, that ty(x,t) obeys the equation 
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^M = _ D ^) + ^W lMl , t) , (4.3) 

which looks like a Schrodinger equation with a harmonic time-dependent potential. 
Eq. (14. 3p is to be solved subject to boundary condition ty(x,T) = 1 for any x. 
We seek the solution of Eq. (14.31) for arbitrary u(t) in the form 

y(x,t) = f(t)ex P (-lg(t)x 2 ) , (4.4) 



2' 



where 



f = Dfg ,/(t = T) = l, (4.5) 



and 



g = 2Dg 2 - w . g(t = T) = 0. (4.6) 

aL> J TUyTjCLT 

Next, we get rid of the nonlinearity in Eq.f l4.6p introducing a new function h obeying 

The function h(t) is solution of the linear second-order differential equation 

h- P& h = 0, (4.8) 
d f Tuj{r)dr 

which has to be solved subject to the boundary conditions 

h{t = T) = 1 , h(t = T) = . (4.9) 
Once h(t) is found, f(t) is determined by f(t) = l/y/Uf) and G(a) by G(a) = f(0) = 

l/y/h(t = 0). 

The moment- generating function for a^2. 

We focus now on the particular case of the weight function u(t) defined by Eq. (11.41) 
with a^2. In this case Eq. (14.81) reads 

(to + t) 
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with 



d Jo r{t + r)~ a dr 
Solution of Eq.f l4.10j) has the form 



> 



(4-11) 



h{t) = Vh + t ]pj v (2vy/a{t Q + t) 2 - a a^j + C 2 K V (2uy/a(t + t) 2 ~ a ff)] , (4.12) 
where and K^(-) are the modified Bessel functions [51], the exponent v is given by 



2- a 



(4.13) 



while the constants C\ and C 2 are chosen to fulfil the boundary conditions in Eqs. fl4.9p . 
so that 



d = 2v^a(t + T) 1 - a o-K u _ 1 (2uy/a{t + T) 2 - a a) 



(4.14) 



and 



C 2 = 2vy/a(t + Ty-^o-I^ (2vy/a(t + T) 2 ~ a o^j . 



(4.15) 



Consequently, we find that h(t = 0) obeys 
h{t = 0) = 



> (o-l)/2 
e \ 7TVZ\/0~ 



1 + e 



2 sin (ttu) 



u-l 



1 + e 



l-a/2 



- I v (vZyfo) I 



l-i/ 



1 + e 



l-a/2 



where e = to/T, and 



(4.16) 



8(1 -a)(2-a)e 2 ~ 



d(e 2 ~ a - (a + e- 1)(1 + C) 1 "") 



(4.17) 



Turning finally to the limit e — >■ 0, we find that the leading small-e behaviour of the 
moment-generating function is given by 
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2 



-d/2 



(4.18) 



for a < 2, while for a > 2 it obeys 



a 



r a - u) I — u 2 

X2j \ VX2 




-d/2 



(4.19) 



where 



d(2-a) , d(a-2) 

Yi = and y 2 = — ; r 

A1 2 A2 2(a-l) 



(4.20) 



The moment- generating function for a = 2. 

We focus next on the particular case a = 2 and consider for convenience a slightly more 
general form of the weight function u(t) by introducing an additional parameter £. We 
stipulate that u(t) is the piece- wise continuous function 



2£/t 2 , for < t < t , 
1/t 2 , fort <t<T. 



^)=rr ; : (^i) 



We seek now the corresponding moment-generating function $(cx) and optimise it in 
what follows considering £ as an optimisation parameter. 

The differential Eq. (14.81) has to be solved now for two intervals t G [0, to] and 
t G [to,r] hi which the "potential" has two different forms. For the first interval, i.e., 
when t G [0,t ], the general solution of Eq. (14. 8[) obeys 

h(t) = ci exp y-^/2a£ > a —J + c 2 exp (^/^^j , (4.22) 

where C\ and c 2 are coefficients to be determined. The parameter a given by Eq. (14. lip 
now reads 

(4.23) 



d(£ + ln(l/e))- 

For the second interval, i.e., when t belongs to [to^], we have 

/ l (^c 3 t(W)/2 + C4t (H/2 ) ( 42 4) 
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where 

S = y/\ + Aaa , (4.25) 

while the coefficients C3 and C4 are to be found from the boundary conditions in 
Eqs. flUD. Tni s yields 

c 3 = ^T-( 1+5 )/ 2 , (4.26) 



28 



and 



5+1 T ( 5 -l)/2 



c 4 = -—77- T ( ]l . (4.27) 
26 

Further on, we require the continuity of h(t) and its first derivative at t = to- We find 
then that 



Cl 



and 



c 2 



(5 + 1) exp {yftE&) ( 5-1 5 5-1 ( 6) \ 
(5+l)e W (-V2E^) f 5-1 s 5-1 , 6) \ 



Consequently, we find that in this case the moment-generating function is given for 
arbitrary e explicitly by 



Ha) = (c 1 + c 2 )- d/2 



(5 + 1) ((i , 5-1 



-e 5 ) cosh {sf2a£o} + 



L 2 ^(<5-i)/2VV 5 + 
5-1 

2^/2aJa 

Now, we are equipped with all necessary results to determine the variance of the 
distribution P( well as the distribution itself. 



A _ 1 / n \ -1 -d/2 



5. The variance of the distribution P{u a ). 

In this section we analyse the behaviour of the variance Var(-u Q ) of the estimator in 
Eq. (11. 31) . First, we calculate exactly the limiting small-e form of Var(ii a ) for arbitrary 
a^2. Further on, we focus on the case a = 2 and determine Var(w Q= 2) for arbitrary e 
and £, Eq. (I4.2ip . We show next that the variance has a minimum for a certain amplitude 
£ = £o P t and present a corresponding expression for the optimised variance. Finally, we 
consider the case when the continuous-space and -time trajectory is recorded at arbitrary 
discrete time moments tj and calculate exactly the weight function u(t) which provides 
the minimal possible variance. 
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5.1. The variance for a^2 and e = 0. 



The variance Vai{u a ) is obtained by differentiating Eqs. (I4.18P or (I4.19P twice with 
respect to o~ and setting a equal to zero. For arbitrary a^2 the variance is then given 
explicitly by 



Var(ri a ) 




2a -3' 



for a < 2, 
for a > 2. 



(5.1) 



The result in the latter equation is depicted in Fig. [2] and shows that, strikingly, the 
variance can be made arbitrarily small in the leading order in e by taking a gradually 
closer to 2, either from above or from below. The slopes at a = 2 + and a = 2~ appear 
to be the same, so that the accuracy of the estimator will be the same for approaching 
a = 2 from above or from below. Equation (15.11) . although formally invalid for a = 2, 
also suggests that the estimator in Eq. (II. 3p with a = 2 has vanishing variance and thus 
possesses an ergodic property. 




Figure 2. (Color online) The variance of the distribution P(u a ) in d = 3 for different 
values of a. Solid line - Eq. (|5.1[) with a < 2 and the dashed line - Eq. (|5.ip with 
a > 2. The symbols correspond to the values obtained from numerical simulations of 
3D random walks for different e, as indicated by the labels. 



A word of caution is now in order. As a matter of fact, we deduce from 
Eq. (I4.16P that finite-e corrections to the result in Eq. (15. ip are of order of 0(e 2 ~ a ) 
for 1 < a < 2, which means that this asymptotical behaviour can be only attained 
when e exp (—1/(2 — a)). In other words, the variance can be made arbitrarily small 
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by choosing a closer to 2, but only at expense of either decreasing the experimental 
resolution time to or increasing the observation time T, which is clearly seen from the 
numerical data shown in Fig. |2j 



5.2. The variance for a = 2 and arbitrary e. 

Differentiating Eq. ( I4.30P with respect to a twice, we find that for arbitrary e the variance 
of the distribution P(u2) is given explicitly by 



£ 31n(l/e)-3(l-e) + 2(l - e)f + £ 2 
Sd~ (£ + m(l/e)) 5 



var(« 2 ) = ^ 1 7 ; : 7:.T:2 ;w . (5.2) 



Notice now that Var(«2) in Eq. ( 15. 21) is a non-monotonic function of £ which attains its 
minimal value when 

c-c (2 + e)ln(l/e)-3(l-e) 

* " e ° pt " ln(l/e) + e-l ' ( j 

This is a somewhat surprising result which states that the optimal choice of the 
amplitude £ in Eq. (I4.2ip . which defines the weight function u(t), actually depends 
on both the time lag to an d on the observation time T. In other words, in order to 
make a proper choice of the amplitude £, one has to know in advance the time through 
which the trajectory is observed. Similar dependence of the optimal parameters on 
the observation time has been recently reported in Refs. [56| [57] , which optimised the 
number of distinct sites visited by intermittent random walks. Note that £ opt — > 2 as 
e — > 0, but for any finite e it is greater than 2. 

Plugging the expression in Eq. (I5.3P into the Eq. ( 15. 2 j) we obtain the corresponding 
optimised variance 

„ ( 31n(l/e) -4 + 5e-e 2 

Var opt (« 2 ) - - ln(1/e) (ln(l/ e ) + l + 2 e )-3(l- e ) ' (5 ' 4) 

From Eq. (ED we find that in 3D Var opt (n 2 ) w 0.144 for e = 10" 3 , Var opt (u 2 ) « 0.096 
for e = 10 -5 and Var opt (M 2 ) ~ 0.082 for e = 10~ 6 . When e — >• 0, Var opt (n 2 ) vanishes in 
a logarithmic way with e at leading order: 

VaroptM-^^. (5.5) 

Therefore, the variance can be made arbitrarily small but at expense of a progressively 
higher resolution or a larger observation time. Note that this is the only case (a = 2) 
in which the estimator defined by Eq. 1 11. 3ft is ergodic. 
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5. 3. Optimisation of the variance for continuous-time and -space trajectories recorded 
at discrete time moments. 

Finally we consider the estimation of the diffusion constant Df when one has a set of 
N temporal points tj such that < to < ^i < ' ' ' < ijv-i = T an d a value Bf , (which 
is one of the components of d- dimensional Brownian motion), for each of these points. 
We consider the least-squares estimator 



N-l 

l 



N-l 
3=0 

where the normalisation is now adsorbed into the weight function Uj, so that 

N-l 

-5>,-E{B*} = l (5.7) 

3=0 

As in the previous subsection, we seek an optimal weight function Uj which minimises 
the variance of the least-squares estimator in Eq. (15. 6ft . Remarkably, this problem can 
be solved exactly for any arbitrary set {tj}. 

The variance of this estimator can be straightforwardly calculated explicitly, for 
arbitrary Uj, to give 

Var(w dis ) = 2 ^ u j u k {tj A t k f , (5.8) 

j,k 

where (tj A t k ) equals the smallest of two numbers tj and th- 
in order to determine the choice of the Uj which minimises the variance of the 
estimator, we minimise the functional 

j,k \ j=0 J 



where A is a Lagrange multiplier enforcing Eq. (15 Jy . Minimising with respect to each 
Uj gives the system of linear equations 

Y^(tj ^hfcoj = Xt k . (5.10) 



To solve this system of equations exactly, we define 

Q fc = ^ Uj , (5.11) 

3>k 
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or, equivalently, 

cjj = flj — fij+i , (5-12) 

for < j < N — 2. Also clearly we have that Qn-i — Wn-i which is compatible with 
defining £7^ = 0. Therefore Eq. f)5.10p can be written as 

- n j+l )t] + t 2 k n k = Xt k . (5.13) 

j<k 

Now subtracting Eq. (15 . 1 3[) for k from the same equation for k + 1 gives the solution 



n k+1 = - — — -, (5.14) 

valid for < k < N — 2, which implies that 

n k = A (5.15) 

which is valid for 1 < A; < iV — 1. In addition, if we set k = in Eq. ( 15. 13ft we find that 

= ^, (5.16) 

which is compatible with Eq. ( 15 . 1 5[) upon defining an additional point t_i = 0. We thus 
find that for 1 < j < N - 2, 

^-^' = (^-o-(^r (5 - 17 » 



while 



and 



^m^t (5 - 18) 



wat„i = Q N _ X = - ^— . (5.19) 

%-l + EjV-2 



The normalisation constraint, Eq. (15. 7p . then yields A as 
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X = ( - + tN ~ 1 + V ^fe+i-^-i) \ / 5 20 n 

\t + h t N -! + t N _ 2 (t j+1 + tfifo + ^_i) y 

Finally the minimal variance can be computed by multiplying Eq.f l5.10j) by uj k and 
summing over k which gives 

(tj A t k ) 2 Uk = A tj — X, (5-21) 

and hence, 

Var(u dis ) = 2 ^ ^ Wfcfo A t k f = 2 A , (5.22) 

j, k 

where we have again used the normalisation condition Eq. (15. 7p . Equations (15. 17)) to 
( 15.20)) define the exact solution of the problem of the optimal estimator for Brownian 
trajectories recorded at discrete time moments. 

In order to compare the results with the continuum case we take the first point to 
to be fixed and write tj = to + A(j — 1) for j > 0, where A is a constant time step, 
A = Tj (N — 1). This gives the following expression for the Lagrange multiplier, which 
defines the variance of the distribution, 

x= ( *° | T h2A y *o + (i-i)A V 1 

l v 2t + A^2T-A + ^(2t + (2j-l)A)(2t + (2j + l)A) y / " ^ > 

Turning to the limit A — > and N —> oo, but keeping the ratio T = N/A fixed, the 
sum in the latter equation becomes a Riemann integral and we find 

,i 1 [ T dt 1, ,T. , 

X-i = l + -^ J = l + - H -), (5.24) 

so that in the leading in e order, for <i-dimensional systems, 

Var(w dis ) = — i— -— . (5.25) 

V J d(2 + ln(l/e)) V ' 

Note that for ln(l/e) ^> 2, the latter equation exhibits exactly the same asymptotic 
behaviour in the limit e — > 0, as the result of the previous Section [5"72| Eq. (15.51) . 

6. Tails of the distribution P(u a ) in d dimensions 

Exact expressions for the moment-generating function allow us to deduce the asymptotic 
behaviour of the distribution P(u a ) for u a < 1 and u a 3> 1. 
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6.1. Asymptotic behaviour of P(u a ) for a^2 and e = 0. 

Large- and small-Mo, asymptotics of P(u a ) can be deduced directly from Eqs. (14.181) 
and (I4.19p . Let us first focus on the small-M Q behaviour, which stems from the large-a 
asymptotical behaviour of the moment-generating function. For a < 2 the latter obeys 



®(a) ~ ^(1+2/(2-00)78 exp _ (6 Vj 




Inverting Eq. ( 16.1 j) we find that for u a — > and a < 2 the distribution P(u a ) shows a 
singular behaviour: 

/ x ( d 2 1 \ 1 , , 

P(u a ) ~ exp -r , (6.2) 



4*1 w 

where the exponent £ is given by 

C=- + -i- eL -r. (6-3) 

s 2 4 12 — Qf| ' V ; 

and the parameter xi is defined in Eq. ( 14.20p . The analogous left tail for a > 2 case 
can be obtained from Eq. ( 16. 2 p by simply making the replacement xi ~+ X2- 

Note that Eq. (I6.2p describes a bell-shaped function whose maximum gives an 
approximation to the most probable value of the estimator u 

2d , lS 

ad + b\2 — a\ 

Note that for any fixed d and a — > 2, the most probable u* — > 1, i.e. to the ensemble 
average value of the estimator in Eq. (jl.3p . Therefore, the least-squares estimators 
outperform the naive end-to-end estimator of the diffusion coefficient, whose distribution 
is given in Eq. (jl.2p and has a bell-shaped form only for d > 3. 

Next, we turn to the large-w Q asymptotical behaviour of the distribution function. 
To do this, it is convenient to rewrite the result in Eq. (14. 18j) as 



= II ( X + 8a ^ 2 - a )*ti-l J "" /2 ' ( 6 - 5 ) 

771=1 



where 7^ jm is the m-th zero of the Bessel function J^z), organised in an ascending order 
[SI]. The large-Ma behaviour of the distribution function stems from the behaviour of 
the moment-generating function in the vicinity of the singular point on the negative 
cr-axis which is closest to a = (all singularities are all located on the negative cx-axis). 
This yields, for a < 2, an exponential decay of the form 
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In a similar fashion, we get that for a > 2 the moment-generating function can be 
represented as 



= II i 1 + 8 (« - 1 W(« - 2 )^,m) 7 , (6-7) 

m=l 

so that in this case the right tail of P(u a ) follows 

P(u a ) ~ uf - 1 exp (-^i • u a ) • (6.8) 

To summarise the results of this subsection, we note the following: 

• when a — > 2, either from above or from below, the small-u a behaviour of P(u a ) 
becomes progressively more singular and small values of u a become very unlikely 
since both xi an d Xi tend to zero and the exponent ( diverges. 

• when a —j- 2, either from above or from below, the inverse relaxation "lengths" 
(i.e., the terms (a — 2) r y 2 _ ul and (2 — a)7 2 „ 1 1 in the exponentials in Eqs. (16. 6p and 
( 16.81) ) diverge, suppressing large values of u a in the distribution. 

Since P(u a ) is normalised for arbitrary a, so that the area below the curve is fixed, 
this implies that P(u a ) tends to the delta-function as a — > 2. 

6.2. Asymptotic behaviour of P(u a ) for a = 2 and small e. 

We focus first on the left tails of the distribution for a = 2 and fixed small e. From 
Eq. (14.301) we find that in the limit a -> oo (so that 5 in Eq. (I4.25P is 5 > 1), fixed 
sufficiently small e, the moment-generating function obeys 




from which equation we readily obtain the following singular form: 

P{U2) ^ exp (-^m.i) i , (6 . 10) 



U2j nil 1 

which holds for u 2 1. Within the opposite limit, i.e., for u 2 ^> 1, the leading 
behaviour of the distribution P(u2) is dominated by the closest to the origin root of the 
denominator in Eq. (I4.30p . Some algebra gives that for e — > the distribution function 
P{u-i) has the following simple form 
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P(u 2 ) ~V exp ^ ^— w 2 J , (6.11) 

where £ opt is the optimised amplitude in Eq. (15 M and xq, in the limit e — > 0, is the root 
of the equation 



_ 2x|\ 1/2 Xq COS (xq) 



{opt/ sin(xo) 2 + ln(l/e) 
The asymptotic behaviour of Xq can be readily obtained: 



(6-12) 





zo = \/^r|i / , K + o ( - 4 ; - - ) 1 . (6.13) 

fopt cot 2 (v^W2j (2 + ln(l/e)) 2 

Therefore, the characteristic decay lengths of the distribution Pin?) from both sides 
from the maximum vanish as l/ln(l/e) when e — > 0. 

7. The distribution P{u a ) in d dimensions 

We turn now to the inversion of the Laplace transform in Eq. (13.11) in order to visualise 
the full distribution P{u a ) and to get an idea of the location of most probable values of 
the estimators in Eq. (11.31) . 

7.1. Inversion of the Laplace transform for a/2 

As we have already remarked, all poles of the moment-generating function $(cr) lie 
on the complex plane on the negative real cx-axis, as can be readily seen from the 
representations in Eqs. (16. 5 j) and (16 .7|) . Setting then 7 = in Eq. (13. 5p . we find 

1 [°° dz cos (zu a — d<p a (z)/2) 



7T 



p(ua) = - 1 v ; /4 , / av " \ (7.1) 



pd (z) 



where, for a < 2, 

Pa(z) = T 2 („) (^f ' (ber 2 ^ ( 2 ^ + bei 2 ^ (2 ^)) , (7.2) 
and the phase is given by 

= arctg (ber^ ^2 J^j /ber„_i ^2 J^A ) , (7.3) 



while for a > 2 we have 



X2\- u L 2 /„ / z \ , , . 2 /„ / ^ 



p Q (z) = r 2 (l-^) ^ beri J 2 . - + beii A 2 , - , (7.4) 




z J \ V V XiJ V V X2 

and 



21 



Distribution of the least-squares estimators 



D. Boyer et al 



(f> a (z) = arctg (ber_ u (2 J^j /ber_„ (^J^j J , (7-5) 

where ber^(x) and bei M (x) are the Kelvin functions [51]. Equation (17. ip defines the 
probability distributions P(u a ) in the leading in e order for arbitrary and arbitrary 

spatial dimension d. 

• — > — ' — 1 — ' — ' — ' — ' — 1 — ' — 1 — 1 — ' — 1 — ' — ' — 1 — ' — 1 — 1 — ' — 1 — 1 — 

) \ 



0.5 1 1.5 2 2.5 




U a U a 

Figure 3. (Color online) The distribution P(u a ) in 3D systems. Upper panel: Colour 
density map of P(u a ) as a function of a, obtained from numerical simulations of 3D 
random walks, with e = 10~ 5 . The solid knots indicate, for different values of a, the 
position of the most probable value of the estimator u a . Lower panel: The distribution 
P{u a ) for different a^2 and with e = 0, obtained by numerical inversion of Eq. (|7.1j) 
for a < 2 (left lower panel) and of Eq. flTTTJ for a > 2 (right lower panel). The symbols 
in the left panel correspond to numerical simulations for (from dark to light), a = — 1 
(circles), a = 1 (squares), a = 3/2 (triangles), and a = 1.95 (stars), and e = 10~ 5 , 
except for a = 1.95 for which we used e = 10~ 7 . 

In Fig. [3] we plot P(u a ) from Eq. (17. ip for a / 2 and e = for three-dimensional 
systems together with the results of numerical simulations. One notices that the 
theoretically predicted distribution P(u a ) becomes more narrow and its maximal value 
grows when a moves towards a = 2, either from above or from below. When a 
approaches 2 from below, the most probable value moves toward the ensemble average 
value (= 1) of the estimator and then starts to move back when a overpasses 2 and grows 
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further. Similarly to the behaviour of the variance, we observe a very good agreement 
between our theoretical predictions and numerical data for a not too close to 2 for 
e = 10 -5 (lower left panel of Figj3]). For a = 1.95 and e as small as 10~ 7 , we still see a 
discrepancy between the numerical data and the asymptotic form of P(u a ) in Eq. ( 17.11) . 
Note that the same slow convergence to zero variance was observed in Figf2]as e — > 0. 




Figure 4. (color online) The distribution P(uj a ) in Eq. (|3.7j) for different a < 2 in 
3D systems. Symbols are the results of numerical simulations. The dashed curve is 
the corresponding result for the end-to-end estimate of the diffusion coefficient in Eq. 



In FigJUwe present the results of numerical simulations for the distribution in Eq. 
(13. 7p of the random variable u defined in Eq. (13. 6p . One notices that as a — > 2, the 
distribution becomes progressively narrower and the peak at u = 1/2 becomes more 
pronounced, which means that it becomes progressively more likely that the diffusion 
coefficients, deduced from two different realisations of Brownian trajectories using the 
weighted least-squares estimators, will have the same value. 



7.2. Two-dimensional systems: series representation of the distributions P(u a ) and 



We proceed further and focus now on the case d = 2 when $(cr) in Eqs. (16. 5p and (16.71) 
have only simple poles. In this case, we readily find via standard residue calculus 



P(u a ) 



2 -u 



oo v 



E 



r(z/ +1) ^ J v (ju-l,m) 
v ' m=l v ' ' 



exp 



(2 - 



-u. 



(7.6) 



for a < 2, while for a > 2 we get 
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P(u r 



E 



i—v,m 



r(2 - V) ^ J^ v (l-v, m ) 
v ' m=l v ' ' 



exp 



(a - 1)4 



(7.7) 



These expressions can be readily plotted using, e.g., Mathematica, and they show exactly 
the same qualitative behaviour (apart from a slightly larger variance), as our results 
obtained by the inversion of the Laplace transform for three-dimensional systems (see, 
Fig. ED. 

We finally present explicit results for the distribution P(ou a ), Eq. (11.41) . for two- 
dimensional systems. For P{u a ) in Eq. (17. 6ft one has, using the definition in Eq. (13. 7ft . 



PM 



2 oo 

E 



E 



7 



v— l.m 



r (l/ + 1) (1 - W a ) J ~_J V (7f-l,n) ^ 

/"°° J / (2-a) / w a 2 2 \ 

x y u Q du a exp I I - _^ j„_ 1>n + j v _ 1>m I u a 



2 00 



(2-a)r(i/ + l)J ^J,(7,-i,n) 
7^ 1 1 



7,- 



l,n 



X 



E 



a J " (T»-i.™) (w«T^-i,„ + (1 - 
Further on, making use of the following equality 



(7.8) 



1 d 



%-l,m duj o 



"V 1 
iv— l.m 



-1 



(7.9) 



and of the definition of the moment-generating function, Eq. ( 14. 181) . 

/"OO 

&(o~) = / P(u a ) exp (— crw Q ) 



(T 



2- a 

2 2-u 



2 



-r -1 



lu-i 2 



(T 



2- a 



'" 1,n ( + 7^-l,n ) » ( 7 - 10 ) 

(2 - a) T(z/ + 1) ^ J v (7v_i,„) V 2 - a 

we can perform one of the summations in Eq. ( 17.81) and, finally, expressing the Bessel 
functions in terms of hypergeometric series, we find 



%- 



l.n 



4a 



+ 7 



d 



^> = ^E 



m=l 



7iLl,m 0F1 [V + 1 



72-1, 



3 Fi 1/. 



1 - CJ a 7„-i, 



Co',. 



(7.11) 



In a similar fashion, for the case a > 2 we obtain 
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4(a-l) d 
a — 2 du r 



m=l 



H v , m \ ( 1 -uj a 7 

7-„,m 0^1 ( 2 - i/, ^- J o^i I 1 - ^ 



q t—v'm 



•(7 



One may readily notice that these two latter distributions are normalised. 

1.95 




Figure 5. (color online) The distribution P(u] a ) in Eq. (|7.11[) for different a < 2 in 
2D systems. 



In Fig. [5] we plot the distribution in Eq. 17.111 for different values of a < 2. 
One notices that similarly to the 3D case, as a — > 2 the distribution becomes 
progressively narrower and the peak at u = 1/2 becomes more pronounced, which 
means that it becomes progressively more likely that the diffusion coefficients, deduced 
from two different realisations of Brownian trajectories using the weighted least-squares 
estimators, will have the same value. One notices as well that for the unweighted LSE 
(a = 1), the distribution P(cj) is rather flat, pointing at large discrepancies between the 
estimates obtained from different trajectories. 



8. Conclusions 



To conclude, we have studied the distribution function P(u a ) of the estimators u a ~ 
T -1 J Q T ui(t) dt, which optimise the least-squares fitting of the diffusion coefficient Df 
of a single <i-dimensional Brownian trajectory B t . We pursued here the optimisation 
further by considering a family of weight functions of the form u(t) = (to + t)~ a , where 
to is a time lag and a is an arbitrary real number, and seeking such values of a for 
which the estimators most efficiently filter out the fluctuations. We calculated P(u a ) 
exactly for arbitrary a and for arbitrary spatial dimension d, and showed that only for 
a = 2 the distribution P(u a ) converges, as e = t /T — > 0, to the Dirac delta-function 
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centered at the ensemble average value of the estimator. This allowed us to conclude 
that only the estimators with a = 2 possess an ergodic property, so that the ensemble 
averaged diffusion coefficient can be obtained with any necessary precision from a single 
trajectory data, but at the expense of a progressively higher experimental resolution. 
For any a^2 the distribution attains, as e — > 0, a certain limiting form with a finite 
variance, which signifies that such estimators are not ergodic. 
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